Loading key packages
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(sf))
suppressPackageStartupMessages(library(zoo))
library(ggplot2)
library(ggspatial)
suppressPackageStartupMessages(library(spatstat))
library(stars)
# Loading required package: abindNiger_price <- read.csv("C:/Users/Laptop/Documents/Munster/5. Data to Knowledge/Food Prices SAHEL/Data/Price_Data/with_lat_long/Niger_prices_with_lat_long.csv")
Mali_price <- read.csv("C:/Users/Laptop/Documents/Munster/5. Data to Knowledge/Food Prices SAHEL/Data/Price_Data/with_lat_long/Mali_prices_with_lat_long.csv")
Burkina_price <- read.csv("C:/Users/Laptop/Documents/Munster/5. Data to Knowledge/Food Prices SAHEL/Data/Price_Data/with_lat_long/BurkinaF_prices_with_lat_long.csv")Niger_prices1 <- Niger_price %>%
mutate(X = row_number()) %>%
mutate(price = as.numeric(price)) %>%
mutate(date = ymd(date)) %>% filter(date >= as.Date("2010-01-01") & date <= as.Date("2020-12-31"))
Mali_prices1 <- Mali_price %>%
mutate(X = row_number()) %>%
mutate(price = as.numeric(price)) %>%
mutate(date = dmy(date)) %>% filter(date >= as.Date("2010-01-01") & date <= as.Date("2020-12-31"))
Burkina_prices1 <- Burkina_price %>%
mutate(X = row_number()) %>%
mutate(price = as.numeric(price)) %>%
mutate(date = ymd(date)) %>% filter(date >= as.Date("2010-01-01") & date <= as.Date("2020-12-31")) %>% arrange(desc(date))# variables for filtering food types:
food <- c("Maize - Retail", "Millet - Retail", "Sorghum - Retail")
food1 <- c("Maize (white) - Retail", "Millet - Retail", "Sorghum (white) - Retail")
Niger_prices4 <- Niger_prices3 %>% filter(cmname %in% food)
Mali_prices4 <- Mali_prices3 %>% filter(cmname %in% food)
Burkina_prices4 <- Burkina_prices3 %>% filter(cmname %in% food1)remove_nonunique_records <- function(data_frame) {
data_frame2 <- data_frame %>%
group_by(quarter, mktname, cmname)
data_frame3 <- data_frame2 %>% mutate(distinct_prices = n_distinct(price))
data_frame4 <- data_frame3[data_frame3$distinct_prices > 2, ] %>% ungroup()
{
return(data_frame4)
}
}
Niger_prices5 <- remove_nonunique_records(Niger_prices4)
Mali_prices5 <- remove_nonunique_records(Mali_prices4)
Burkina_prices5 <- remove_nonunique_records(Burkina_prices4)calculate_standard_dev <- function(data_frame) {
data_frame2 <- data_frame %>%
group_by(cmname) %>%
mutate(price_norm = (price - min(price)) / (max(price) - min(price))) %>%
ungroup()
data_frame3 <- data_frame2 %>%
group_by(quarter, mktname, cmname) %>%
mutate(stdev = sd(price_norm, na.rm =TRUE)) %>%
ungroup()
{
return(data_frame3)
}
}
Niger_prices6 <- calculate_standard_dev(Niger_prices5)
Mali_prices6 <- calculate_standard_dev(Mali_prices5)
Burkina_prices6 <- calculate_standard_dev(Burkina_prices5)calculate_mean_stdev <- function(data_frame) {
data_frame2 <- data_frame %>%
group_by(quarter, mktname) %>%
mutate(stdev_mean = mean(stdev, na.rm = TRUE)) %>%
ungroup()
{
return(data_frame2)
}
}
Niger_prices7 <- calculate_mean_stdev(Niger_prices6)
Mali_prices7 <- calculate_mean_stdev(Mali_prices6)
Burkina_prices7 <- calculate_mean_stdev(Burkina_prices6)merge_same_location_markets <- function(data_frame) {
data_frame2 <- data_frame %>%
group_by(quarter, lat, long) %>%
mutate(stdev_mean = mean(stdev_mean, na.rm =TRUE)) %>%
ungroup()
{
return(data_frame2)
}
}
Niger_price8 <- Niger_prices7 #does not have any location markets
Mali_price8 <- merge_same_location_markets(Mali_prices7)
Burkina_price8 <- Burkina_prices7 #does not have any location marketsmean_diff_column <- function(data_frame) {
data_frame2<- data_frame %>%
group_by(quarter, cmname) %>%
mutate(avg_nat_mean= mean(price)) %>%
ungroup()
data_frame3<- data_frame2%>%
group_by(quarter, cmname, mktname) %>%
mutate(avg_local_price= mean(price)) %>%
ungroup()
data_frame4 <-data_frame3 %>% mutate(Percent_Difference = (avg_local_price - avg_nat_mean)/avg_nat_mean*100)
data_frame5<-data_frame4%>%
group_by(quarter, mktname) %>%
mutate(avg_difference = mean(Percent_Difference)) %>%
ungroup()
{
return(data_frame5)
}
}
Niger_price10 <- mean_diff_column(Niger_price8)
Mali_price10 <- mean_diff_column(Mali_price8)
Burkina_price10 <- mean_diff_column(Burkina_price8)
# add a "quarter" column
Niger_conflicts$quarter = Niger_conflicts$event_date
Mali_conflicts$quarter = Mali_conflicts$event_date
Burkina_conflicts$quarter = Burkina_conflicts$event_date
# function for pre-processing the conflict data
conflicts_filter = c("Battles", "Violence against civilians", "Explosions/Remote violence")
#conflicts_filter = c("Battles", "Protests", "Violence against civilians", "Explosions/Remote violence", "Riots")
Pre_process_conflicts <- function(data_frame) {
data_frame2 <- data_frame %>% mutate(event_date = dmy(event_date)) %>%
rename("date" = "event_date", "type" = "event_type", "lat" = "latitude", "long" = "longitude") %>% filter(type %in% conflicts_filter)
data_frame2$date <- as.yearmon(data_frame2$date, "%Y%m")
data_frame2$quarter <- as.yearqtr(dmy(data_frame2$quarter), "%Y%m")
final_pre_processsed <- data_frame2 %>% arrange(date)
{
return(final_pre_processsed)
}
}
# pre-process the conflict data
Niger_conflicts3 <- Pre_process_conflicts(Niger_conflicts)
Mali_conflicts3 <- Pre_process_conflicts(Mali_conflicts)
Burkina_conflicts3 <- Pre_process_conflicts(Burkina_conflicts)WKT_conversion_function <- function(data_frame) {
data_frame$geometry <- paste("POINT(",data_frame$long, ",", data_frame$lat,")")
data_frame2 <- data_frame %>% mutate(geometry = str_replace_all(geometry, " ", ""))
data_frame3 <- data_frame2 %>% mutate(geometry = str_replace_all(geometry, ",", " "))
{
return(data_frame3)
}
}
WKTcol_to_sfc_and_make_sf_collection <- function(data_frame) {
data_frame$geometry <- st_as_sfc(data_frame$geometry)
final_sf_collection <- st_sf(data_frame, crs = "EPSG:4326")
{
return(final_sf_collection)
}
}
Niger_conflicts4 <- WKT_conversion_function(Niger_conflicts3)
Niger_conflicts5 <- WKTcol_to_sfc_and_make_sf_collection(Niger_conflicts4)
Mali_conflicts4 <- WKT_conversion_function(Mali_conflicts3)
Mali_conflicts5 <- WKTcol_to_sfc_and_make_sf_collection(Mali_conflicts4)
Burkina_conflicts4 <- WKT_conversion_function(Burkina_conflicts3)
Burkina_conflicts5 <- WKTcol_to_sfc_and_make_sf_collection(Burkina_conflicts4)
Niger_price11 <- WKT_conversion_function(Niger_price10)
Niger_price12 <- WKTcol_to_sfc_and_make_sf_collection(Niger_price11)
Mali_price11 <- WKT_conversion_function(Mali_price10)
Mali_price12 <- WKTcol_to_sfc_and_make_sf_collection(Mali_price11)
Burkina_price11 <- WKT_conversion_function(Burkina_price10)
Burkina_price12 <- WKTcol_to_sfc_and_make_sf_collection(Burkina_price11)Niger_country <- st_read("C:/Users/Laptop/Documents/Munster/5. Data to Knowledge/Food Prices SAHEL/Data/Country_Boundaries/Niger_boundaries.geojsonl.json")
# Reading layer `Niger_boundaries.geojsonl' from data source
# `C:\Users\Laptop\Documents\Munster\5. Data to Knowledge\Food Prices SAHEL\Data\Country_Boundaries\Niger_boundaries.geojsonl.json'
# using driver `GeoJSON'
# Simple feature collection with 1 feature and 2 fields
# Geometry type: MULTIPOLYGON
# Dimension: XY
# Bounding box: xmin: 0.1529411 ymin: 11.69577 xmax: 15.97032 ymax: 23.51735
# Geodetic CRS: WGS 84
Mali_country <- st_read("/Users/Laptop/Documents/Munster/5. Data to Knowledge/Food Prices SAHEL/Data/Country_Boundaries/Mali_boundaries.geojsonl.json")
# Reading layer `Mali_boundaries.geojsonl' from data source
# `C:\Users\Laptop\Documents\Munster\5. Data to Knowledge\Food Prices SAHEL\Data\Country_Boundaries\Mali_boundaries.geojsonl.json'
# using driver `GeoJSON'
# Simple feature collection with 1 feature and 2 fields
# Geometry type: MULTIPOLYGON
# Dimension: XY
# Bounding box: xmin: -12.26413 ymin: 10.14005 xmax: 4.235638 ymax: 24.99506
# Geodetic CRS: WGS 84
Burkina_country <- st_read("C:/Users/Laptop/Documents/Munster/5. Data to Knowledge/Food Prices SAHEL/Data/Country_Boundaries/Burkina_boundaries.geojsonl.json")
# Reading layer `Burkina_boundaries.geojsonl' from data source
# `C:\Users\Laptop\Documents\Munster\5. Data to Knowledge\Food Prices SAHEL\Data\Country_Boundaries\Burkina_boundaries.geojsonl.json'
# using driver `GeoJSON'
# Simple feature collection with 1 feature and 2 fields
# Geometry type: MULTIPOLYGON
# Dimension: XY
# Bounding box: xmin: -5.522578 ymin: 9.391883 xmax: 2.390169 ymax: 15.07991
# Geodetic CRS: WGS 84
All_countries <- st_read("C:/Users/Laptop/Documents/Munster/5. Data to Knowledge/Food Prices SAHEL/Data/Country_Boundaries/All_countries_boundaries.geojsonl.json")
# Reading layer `All_countries_boundaries.geojsonl' from data source
# `C:\Users\Laptop\Documents\Munster\5. Data to Knowledge\Food Prices SAHEL\Data\Country_Boundaries\All_countries_boundaries.geojsonl.json'
# using driver `GeoJSONSeq'
# Simple feature collection with 3 features and 2 fields
# Geometry type: MULTIPOLYGON
# Dimension: XY
# Bounding box: xmin: -12.26413 ymin: 9.391883 xmax: 15.97032 ymax: 24.99506
# Geodetic CRS: WGS 84
ggplot() +
geom_sf(data = All_countries, fill= "antiquewhite") +
annotation_scale(location = "br", width_hint = 0.5) +
annotation_north_arrow(location = "tl", which_north = "true", pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"), style = north_arrow_fancy_orienteering) +
xlab("Longitude") +
ylab("Latitude") +
ggtitle("Countries of Interest") +
theme(plot.title = element_text(size=20, face="bold")) +
theme(panel.grid.major = element_line(color = gray(.5), linetype = "dashed", size = 0.5), panel.background = element_rect(fill = "aliceblue")) +
annotate(geom = "text", x = -2, y = 18, label = "Mali", fontface = "bold", color = "grey22", size = 4) +
annotate(geom = "text", x = 10, y = 17, label = "Niger", fontface = "bold",color = "grey22", size = 4) +
annotate(geom = "text", x = -1, y = 12, label = "Burkina Faso", fontface = "bold", color = "grey22", size = 4)
# Scale on map varies by more than 10%, scale bar may be inaccurateggplot() +
geom_sf(data = All_countries, fill= "antiquewhite") +
geom_sf(data = rbind(Niger_conflicts5, Mali_conflicts5, Burkina_conflicts5), size=0.3, aes(col = "red")) +
annotation_scale(location = "br", width_hint = 0.5) +
annotation_north_arrow(location = "tl", which_north = "true", pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"), style = north_arrow_fancy_orienteering) +
xlab("Longitude") +
ylab("Latitude") +
ggtitle("Spatial Distribution of Conflicts (2010-2020)") +
theme(plot.title = element_text(size=20, face="bold")) +
theme(panel.grid.major = element_line(color = gray(.5), linetype = "dashed", size = 0.5), panel.background = element_rect(fill = "aliceblue")) +
scale_color_manual(values=c("red"), labels=c("Conflicts")) + labs(color='Legend')
# Scale on map varies by more than 10%, scale bar may be inaccurateggplot() +
geom_sf(data = All_countries, fill= "antiquewhite") +
geom_sf(data = rbind(Niger_price12, Mali_price12, Burkina_price12), size=1, aes(col = "navyblue"))+
annotation_scale(location = "br", width_hint = 0.5) +
annotation_north_arrow(location = "tl", which_north = "true", pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"), style = north_arrow_fancy_orienteering) +
xlab("Longitude") +
ylab("Latitude") +
ggtitle("Spatial Distribution of Food Markets (2010-2020)") +
theme(plot.title = element_text(size=20, face="bold")) +
theme(panel.grid.major = element_line(color = gray(.5), linetype = "dashed", size = 0.5), panel.background = element_rect(fill = "aliceblue")) +
scale_color_manual(values=c("navyblue"), labels=c("Markets")) + labs(color='Legend')
# Scale on map varies by more than 10%, scale bar may be inaccurateNiger_mean_quarter <- Niger_price12 %>% group_by(quarter) %>% summarise(stdev_mean = mean(stdev_mean))
Mali_mean_quarter <- Mali_price12 %>% group_by(quarter) %>% summarise(stdev_mean = mean(stdev_mean))
Burkina_mean_quarter <- Burkina_price12 %>% group_by(quarter) %>% summarise(stdev_mean = mean(stdev_mean))
library("plotly")
#
# Attaching package: 'plotly'
# The following object is masked from 'package:ggplot2':
#
# last_plot
# The following object is masked from 'package:stats':
#
# filter
# The following object is masked from 'package:graphics':
#
# layout
ggplotly(
ggplot() +
geom_line(data= Niger_mean_quarter, aes(y=stdev_mean, x=quarter, col="Niger")) +
geom_line(data=Mali_mean_quarter, aes(y = stdev_mean, x=quarter,col = "Mali")) +
geom_line(data= Burkina_mean_quarter, aes(y = stdev_mean, x=quarter,col="Burkina Faso")) +
labs(title="Quarterly Price Standard Deviation per Country", y="standard deviation of min-max normalised prices", x="quarter",
color=NULL) + scale_color_manual(labels = c("Niger", "Burkina Faso", "Mali"),
values = c("Niger"="#00ba38", "Burkina Faso"="#f8766d", "Mali" = "blue")) +
theme(axis.text.x = element_text( vjust=0.5 , size = 8), plot.title = element_text(size=16, face="bold"))
)Niger_quarterly_conflicts2 <- Niger_conflicts3 %>% group_by(quarter) %>% summarise(n_conflicts = n_distinct(data_id))
Mali_quarterly_conflicts2 <- Mali_conflicts3 %>% group_by(quarter) %>% summarise(n_conflicts = n_distinct(data_id))
Burkina_quarterly_conflicts2 <- Burkina_conflicts3 %>% group_by(quarter) %>% summarise(n_conflicts = n_distinct(data_id))
ggplotly(
ggplot() +
geom_line(data = Niger_quarterly_conflicts2, aes(x = quarter, y = n_conflicts, col="Niger")) +
geom_line(data = Burkina_quarterly_conflicts2, aes(x = quarter, y = n_conflicts, col="Burkina Faso")) +
geom_line(data = Mali_quarterly_conflicts2, aes(x = quarter, y = n_conflicts, col="Mali")) +
labs(title="Quarterly Conflicts per Country", y="number of conflicts", x="quarter", color=NULL) +
scale_color_manual(labels = c("Niger", "Burkina Faso", "Mali"),
values = c("Niger"="#00ba38", "Burkina Faso"="#f8766d", "Mali" = "blue")) +
theme(axis.text.x = element_text( vjust=0.5 , size = 8), plot.title = element_text(size=16, face="bold"))
)FINAL_results <- data.frame(mktname = 0, final_std = 0, v = 0, final_avg_diff = 0, geometry = 0)
merged_list_margets <- rbind(Niger_price12, Mali_price12, Burkina_price12)
merged_list_margets2 <- merged_list_margets %>% filter(quarter != "2021 Q1")
unique_quarters <- c(unique(merged_list_margets2$quarter))
# All_countries polygons
All_countries2 <- All_countries %>% st_transform(3857) %>% select(ADMIN)
# Bounding Box for All_countries3
bb <- st_geometry(All_countries2) %>% st_union()
BW_df <- data.frame()
# list of bandwidths for the point density maps
sigma_list <- list(200000, 150000, 100000, 90000, 80000, 70000, 60000 ,50000, 45000, 40000 ,35000, 30000, 25000, 20000, 15000, 10000, 5000)
#
# dataframes for final correlation results per bandwidth
sd_correlations <- data.frame(bandwidths = c(200000, 150000, 100000, 90000, 80000, 70000, 60000 ,50000, 45000, 40000 ,35000, 30000, 25000, 20000, 15000, 10000,5000), correlation = NA)
diff_from_mean_correlations <- data.frame(bandwidths = c(200000, 150000, 100000, 90000, 80000, 70000, 60000 ,50000, 45000, 40000 ,35000, 30000, 25000, 20000, 15000, 10000, 5000), correlation = NA)
#This loop runs all the steps of the analysis one time for each of the bandwidths at a time. For each bandwidth the analysis is also run one time for each quarter between 2010-2020
for (bandwidht in sigma_list) {
for (Q in unique_quarters) {
# Step 1: Kernel Point Density Surfaces
# Point data for conflicts for each country
Niger_points <- Niger_conflicts5 %>% st_transform(3857) %>% filter(quarter == Q)
Mali_points <- Mali_conflicts5 %>% st_transform(3857) %>% filter(quarter == Q)
Burkina_points <- Burkina_conflicts5 %>% st_transform(3857) %>% filter(quarter == Q)
# Format as point pattern data
point_patern = c(bb, st_geometry(Niger_points), st_geometry(Mali_points), st_geometry(Burkina_points)) %>% as.ppp()
skip_to_next <- FALSE
# The density function takes the point pattern data and the bandwidth as arguments. The dimyx refers the pixel resolution
tryCatch(density_map <- density(point_patern, sigma = bandwidht, dimyx = 1000, verbose = FALSE), error = function(e) { skip_to_next <<- TRUE})
if(skip_to_next) { next }
density_map_stars <- st_as_stars(density_map) %>% st_set_crs(3857)
Niger_price13 <- Niger_price12 %>% filter(quarter == Q)
Mali_price13 <- Mali_price12 %>% filter(quarter == Q)
Burkina_price13 <- Burkina_price12 %>% filter(quarter == Q)
# Market points
All_markets_final <- rbind(Niger_price13, Mali_price13, Burkina_price13) %>% st_transform(3857)
# count column
All_markets_final$counts <- 1
# Aggregate std of all commodities to one row per mktname/quarter:
All_markets_final2 <- All_markets_final %>% group_by(mktname) %>% summarise(final_std = sum(stdev_mean/sum(counts)), final_avg_diff = sum(avg_difference/sum(counts)))
# Step 2: Extract values for each market point from density map:
extracted <- st_extract(density_map_stars, All_markets_final2)
Final_join <- st_join(All_markets_final2, extracted, join = st_equals)
FINAL_results <- rbind(FINAL_results, Final_join)
}
# Step 3: Regression analysis. For each bandwidth a correlation coefficient is calculated between the two metrics and teh conflict density extracted at the markets within each of the quarters.
sd_correlations[sd_correlations$bandwidths == bandwidht, "correlation"] <- with(na.omit(FINAL_results), cor(final_std, v))
diff_from_mean_correlations[sd_correlations$bandwidths == bandwidht, "correlation"] <- with(na.omit(FINAL_results), cor(final_avg_diff, v))
# The data frame is reset before running next bandwidth
FINAL_results <- data.frame(mktname = 0, final_std = 0, v = 0, final_avg_diff = 0, geometry = 0)
}PS: Some of the conflict points fall outside the country boundaries used as a bounding box and are therefore not included in the point density estimations. Warnings about duplicated points are ignored as this just means some conflict events have been recorded at the same location. All the points are still included in the density estimation.
Niger_price16 <- Niger_price12 %>% filter(quarter == "2016 Q1")
Mali_price16 <- Mali_price12 %>% filter(quarter == "2016 Q1")
Burkina_price16 <- Burkina_price12 %>% filter(quarter == "2016 Q1")
All_markets_final <- rbind(Niger_price16, Mali_price16, Burkina_price16) %>% st_transform(3857)
plot_bandwidths <- list(200000, 45000, 5000)
#
for (plot_bandwidth in plot_bandwidths) {
# Point data for conflicts for each country
Niger_points <- Niger_conflicts5 %>% st_transform(3857) %>% filter(quarter == "2016 Q1")
Mali_points <- Mali_conflicts5 %>% st_transform(3857) %>% filter(quarter == "2016 Q1")
Burkina_points <- Burkina_conflicts5 %>% st_transform(3857) %>% filter(quarter == "2016 Q1")
# Format as point pattern data
point_patern = c(bb, st_geometry(Niger_points), st_geometry(Mali_points), st_geometry(Burkina_points)) %>% as.ppp()
plot_density_map <- density(point_patern, sigma = plot_bandwidth, dimyx = 1000)
plot(plot_density_map, main= paste("Point density map for Q1 2017 with bandwidth", plot_bandwidth), reset = FALSE)
plot(All_markets_final, size=1, col = "green", add = TRUE)
}ggplotly(
ggplot(sd_correlations, aes(x=bandwidths, y= correlation)) +
geom_smooth()+
geom_point() +
ggtitle("Price Standard Deviation vs. Conflict Events")+
ylab("Correlation coefficients") +
xlab("Bandwidhts (meters)") +
theme(plot.title = element_text(size=16, face="bold"))
)